//
// File: freqz.cpp
//
// MATLAB Coder version            : 5.4
// C/C++ source code generated on  : 21-Mar-2025 17:03:22
//

// Include Files
#include "freqz.h"
#include "FFTImplementationCallback.h"
#include "calc_TimeFreq_Ar_data.h"
#include "rt_nonfinite.h"
#include "coder_array.h"
#include <cmath>

// Function Definitions
//
// Arguments    : coder::array<creal_T, 1U> &in1
//                const coder::array<creal_T, 2U> &in2
//                const coder::array<creal_T, 2U> &in3
// Return Type  : void
//
void binary_expand_op(coder::array<creal_T, 1U> &in1,
                      const coder::array<creal_T, 2U> &in2,
                      const coder::array<creal_T, 2U> &in3)
{
  int i;
  int loop_ub;
  int stride_0_0;
  int stride_1_0;
  if (in3.size(1) == 1) {
    i = in2.size(1);
  } else {
    i = in3.size(1);
  }
  in1.set_size(i);
  stride_0_0 = (in2.size(1) != 1);
  stride_1_0 = (in3.size(1) != 1);
  if (in3.size(1) == 1) {
    loop_ub = in2.size(1);
  } else {
    loop_ub = in3.size(1);
  }
  for (i = 0; i < loop_ub; i++) {
    double ai;
    double ar;
    double bi;
    double br;
    double im;
    double re;
    int br_tmp;
    ar = in2[i * stride_0_0].re;
    ai = in2[i * stride_0_0].im;
    br_tmp = i * stride_1_0;
    br = in3[br_tmp].re;
    bi = in3[br_tmp].im;
    if (bi == 0.0) {
      if (ai == 0.0) {
        re = ar / br;
        im = 0.0;
      } else if (ar == 0.0) {
        re = 0.0;
        im = ai / br;
      } else {
        re = ar / br;
        im = ai / br;
      }
    } else if (br == 0.0) {
      if (ar == 0.0) {
        re = ai / bi;
        im = 0.0;
      } else if (ai == 0.0) {
        re = 0.0;
        im = -(ar / bi);
      } else {
        re = ai / bi;
        im = -(ar / bi);
      }
    } else {
      double brm;
      brm = std::abs(br);
      im = std::abs(bi);
      if (brm > im) {
        double s;
        s = bi / br;
        im = br + s * bi;
        re = (ar + s * ai) / im;
        im = (ai - s * ar) / im;
      } else if (im == brm) {
        double s;
        if (br > 0.0) {
          s = 0.5;
        } else {
          s = -0.5;
        }
        if (bi > 0.0) {
          im = 0.5;
        } else {
          im = -0.5;
        }
        re = (ar * s + ai * im) / brm;
        im = (ai * s - ar * im) / brm;
      } else {
        double s;
        s = br / bi;
        im = bi + s * br;
        re = (s * ar + ai) / im;
        im = (s * ai - ar) / im;
      }
    }
    in1[i].re = re;
    in1[i].im = im;
  }
}

//
// Arguments    : double bin
//                double npoints
//                ::coder::array<creal_T, 1U> &h
// Return Type  : void
//
namespace coder {
void firfreqz(double bin, double npoints, ::coder::array<creal_T, 1U> &h)
{
  array<creal_T, 1U> b_fv;
  array<creal_T, 1U> fv;
  array<creal_T, 1U> wwc;
  array<creal_T, 1U> yCol;
  array<double, 2U> costab;
  array<double, 2U> sintab;
  array<double, 2U> sintabinv;
  int iDelta2;
  int ihi;
  if (!(npoints < 1.0)) {
    int len_tmp_tmp;
    int nt_re_tmp;
    boolean_T useRadix2;
    len_tmp_tmp = static_cast<int>(npoints);
    useRadix2 =
        ((static_cast<int>(npoints) & (static_cast<int>(npoints) - 1)) == 0);
    internal::FFTImplementationCallback::get_algo_sizes(
        static_cast<int>(npoints), useRadix2, &ihi, &iDelta2);
    internal::FFTImplementationCallback::generate_twiddle_tables(
        iDelta2, useRadix2, costab, sintab, sintabinv);
    if (useRadix2) {
      yCol.set_size(len_tmp_tmp);
      if (static_cast<int>(npoints) > 1) {
        yCol.set_size(len_tmp_tmp);
        for (nt_re_tmp = 0; nt_re_tmp < len_tmp_tmp; nt_re_tmp++) {
          yCol[nt_re_tmp].re = 0.0;
          yCol[nt_re_tmp].im = 0.0;
        }
      }
      if (static_cast<int>(npoints) != 1) {
        internal::FFTImplementationCallback::doHalfLengthRadix2(
            bin, yCol, static_cast<int>(npoints), costab, sintab);
      } else {
        double b_nt_re_tmp;
        double nt_im;
        double nt_re;
        double twid_im;
        int i;
        int iDelta;
        int k;
        int nInt2m1;
        int nRowsD2;
        iDelta = static_cast<int>(npoints) - 2;
        nRowsD2 = static_cast<int>(npoints) / 2;
        k = nRowsD2 / 2;
        yCol[0].re = bin;
        yCol[0].im = 0.0;
        if (static_cast<int>(npoints) > 1) {
          for (i = 0; i <= iDelta; i += 2) {
            b_nt_re_tmp = yCol[i + 1].re;
            nt_im = yCol[i + 1].im;
            twid_im = yCol[i].re;
            nt_re = yCol[i].im;
            yCol[i + 1].re = twid_im - b_nt_re_tmp;
            yCol[i + 1].im = nt_re - nt_im;
            yCol[i].re = twid_im + b_nt_re_tmp;
            yCol[i].im = nt_re + nt_im;
          }
        }
        iDelta = 2;
        iDelta2 = 4;
        nInt2m1 = ((k - 1) << 2) + 1;
        while (k > 0) {
          int rt;
          for (i = 0; i < nInt2m1; i += iDelta2) {
            nt_re_tmp = i + iDelta;
            nt_re = yCol[nt_re_tmp].re;
            nt_im = yCol[nt_re_tmp].im;
            yCol[nt_re_tmp].re = yCol[i].re - nt_re;
            yCol[nt_re_tmp].im = yCol[i].im - nt_im;
            yCol[i].re = yCol[i].re + nt_re;
            yCol[i].im = yCol[i].im + nt_im;
          }
          rt = 1;
          for (int nInt2{k}; nInt2 < nRowsD2; nInt2 += k) {
            double twid_re;
            twid_re = costab[nInt2];
            twid_im = sintab[nInt2];
            i = rt;
            ihi = rt + nInt2m1;
            while (i < ihi) {
              nt_re_tmp = i + iDelta;
              b_nt_re_tmp = yCol[nt_re_tmp].im;
              nt_im = yCol[nt_re_tmp].re;
              nt_re = twid_re * nt_im - twid_im * b_nt_re_tmp;
              nt_im = twid_re * b_nt_re_tmp + twid_im * nt_im;
              yCol[nt_re_tmp].re = yCol[i].re - nt_re;
              yCol[nt_re_tmp].im = yCol[i].im - nt_im;
              yCol[i].re = yCol[i].re + nt_re;
              yCol[i].im = yCol[i].im + nt_im;
              i += iDelta2;
            }
            rt++;
          }
          k /= 2;
          iDelta = iDelta2;
          iDelta2 += iDelta2;
          nInt2m1 -= iDelta;
        }
      }
    } else {
      double nt_im;
      double nt_re;
      int iDelta;
      if ((static_cast<int>(npoints) != 1) &&
          ((static_cast<int>(npoints) & 1) == 0)) {
        int nInt2;
        int nInt2m1;
        int rt;
        iDelta2 = static_cast<int>(npoints) / 2;
        nInt2m1 = (iDelta2 + iDelta2) - 1;
        wwc.set_size(nInt2m1);
        rt = 0;
        wwc[iDelta2 - 1].re = 1.0;
        wwc[iDelta2 - 1].im = 0.0;
        nInt2 = iDelta2 << 1;
        for (int k{0}; k <= iDelta2 - 2; k++) {
          iDelta = ((k + 1) << 1) - 1;
          if (nInt2 - rt <= iDelta) {
            rt += iDelta - nInt2;
          } else {
            rt += iDelta;
          }
          nt_im = -3.1415926535897931 * static_cast<double>(rt) /
                  static_cast<double>(iDelta2);
          if (nt_im == 0.0) {
            nt_re = 1.0;
            nt_im = 0.0;
          } else {
            nt_re = std::cos(nt_im);
            nt_im = std::sin(nt_im);
          }
          nt_re_tmp = (iDelta2 - k) - 2;
          wwc[nt_re_tmp].re = nt_re;
          wwc[nt_re_tmp].im = -nt_im;
        }
        nt_re_tmp = nInt2m1 - 1;
        for (int k{nt_re_tmp}; k >= iDelta2; k--) {
          wwc[k] = wwc[(nInt2m1 - k) - 1];
        }
      } else {
        int nInt2;
        int nInt2m1;
        int rt;
        nInt2m1 = (static_cast<int>(npoints) + static_cast<int>(npoints)) - 1;
        wwc.set_size(nInt2m1);
        rt = 0;
        wwc[static_cast<int>(npoints) - 1].re = 1.0;
        wwc[static_cast<int>(npoints) - 1].im = 0.0;
        nInt2 = static_cast<int>(npoints) << 1;
        for (int k{0}; k <= len_tmp_tmp - 2; k++) {
          iDelta = ((k + 1) << 1) - 1;
          if (nInt2 - rt <= iDelta) {
            rt += iDelta - nInt2;
          } else {
            rt += iDelta;
          }
          nt_im = -3.1415926535897931 * static_cast<double>(rt) /
                  static_cast<double>(static_cast<int>(npoints));
          if (nt_im == 0.0) {
            nt_re = 1.0;
            nt_im = 0.0;
          } else {
            nt_re = std::cos(nt_im);
            nt_im = std::sin(nt_im);
          }
          nt_re_tmp = (static_cast<int>(npoints) - k) - 2;
          wwc[nt_re_tmp].re = nt_re;
          wwc[nt_re_tmp].im = -nt_im;
        }
        nt_re_tmp = nInt2m1 - 1;
        for (int k{nt_re_tmp}; k >= len_tmp_tmp; k--) {
          wwc[k] = wwc[(nInt2m1 - k) - 1];
        }
      }
      yCol.set_size(len_tmp_tmp);
      if (static_cast<int>(npoints) > 1) {
        yCol.set_size(len_tmp_tmp);
        for (nt_re_tmp = 0; nt_re_tmp < len_tmp_tmp; nt_re_tmp++) {
          yCol[nt_re_tmp].re = 0.0;
          yCol[nt_re_tmp].im = 0.0;
        }
      }
      if ((ihi != 1) && ((static_cast<int>(npoints) & 1) == 0)) {
        internal::FFTImplementationCallback::doHalfLengthBluestein(
            bin, yCol, static_cast<int>(npoints), ihi, wwc, costab, sintab,
            costab, sintabinv);
      } else {
        double twid_im;
        double twid_re;
        yCol[0].re = wwc[static_cast<int>(npoints) - 1].re * bin;
        yCol[0].im = wwc[static_cast<int>(npoints) - 1].im * -bin;
        for (int k{2}; k <= len_tmp_tmp; k++) {
          yCol[k - 1].re = 0.0;
          yCol[k - 1].im = 0.0;
        }
        internal::FFTImplementationCallback::r2br_r2dit_trig_impl(
            yCol, ihi, costab, sintab, fv);
        internal::FFTImplementationCallback::r2br_r2dit_trig_impl(
            wwc, ihi, costab, sintab, b_fv);
        b_fv.set_size(fv.size(0));
        iDelta = fv.size(0);
        for (nt_re_tmp = 0; nt_re_tmp < iDelta; nt_re_tmp++) {
          twid_im = fv[nt_re_tmp].re;
          nt_im = b_fv[nt_re_tmp].im;
          nt_re = fv[nt_re_tmp].im;
          twid_re = b_fv[nt_re_tmp].re;
          b_fv[nt_re_tmp].re = twid_im * twid_re - nt_re * nt_im;
          b_fv[nt_re_tmp].im = twid_im * nt_im + nt_re * twid_re;
        }
        internal::FFTImplementationCallback::r2br_r2dit_trig_impl(
            b_fv, ihi, costab, sintabinv, fv);
        if (fv.size(0) > 1) {
          nt_im = 1.0 / static_cast<double>(fv.size(0));
          iDelta = fv.size(0);
          for (nt_re_tmp = 0; nt_re_tmp < iDelta; nt_re_tmp++) {
            fv[nt_re_tmp].re = nt_im * fv[nt_re_tmp].re;
            fv[nt_re_tmp].im = nt_im * fv[nt_re_tmp].im;
          }
        }
        nt_re_tmp = wwc.size(0);
        for (int k{len_tmp_tmp}; k <= nt_re_tmp; k++) {
          nt_im = wwc[k - 1].re;
          nt_re = fv[k - 1].im;
          twid_re = wwc[k - 1].im;
          twid_im = fv[k - 1].re;
          iDelta = k - static_cast<int>(npoints);
          yCol[iDelta].re = nt_im * twid_im + twid_re * nt_re;
          yCol[iDelta].im = nt_im * nt_re - twid_re * twid_im;
        }
      }
    }
    h.set_size(len_tmp_tmp);
    for (nt_re_tmp = 0; nt_re_tmp < len_tmp_tmp; nt_re_tmp++) {
      h[nt_re_tmp] = yCol[nt_re_tmp];
    }
  }
}

} // namespace coder

//
// File trailer for freqz.cpp
//
// [EOF]
//
